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ABSTRACT 

The direct collapse model for the formation of massive seed black holes in the early 
Universe attempts to explain the observed number density of supermassive black holes 
(SMBHs) at z ~ 6 by assuming that they grow from seeds with masses M > 10 4 M 0 
that form by the direct collapse of metal-free gas in atomic cooling halos in which 
H 2 cooling is suppressed by a strong extragalactic radiation field. The viability of 
this model depends on the strength of the radiation field required to suppress H 2 
cooling, J cr i t : if this is too large, then too few seeds will form to explain the observed 
number density of SMBHs. In order to determine J cr ; t reliably, we need to be able 
to accurately model the formation and destruction of H 2 in gas illuminated by an 
extremely strong radiation Held. In this paper, we use a reaction-based reduction 
technique to analyze the chemistry of H 2 in these conditions, allowing us to identify 
the key chemical reactions that are responsible for determining the value of J cr it- 
We construct a reduced network of 26 reactions that allows us to determine J cr i t 
accurately, and compare it with previous treatments in the literature. We show that 
previous studies have often omitted one or more important chemical reactions, and 
that these omissions introduce an uncertainty of up to a factor of three into previous 
determinations of J cr it- 

Key words: astrochemistry - hydrodynamics - methods: numerical - molecular 
processes - cosmology: theory - quasars: general 


1 INTRODUCTION 

In recent years, infrared sky surveys such as UKIDSS have 
discovered a number of quasars at redshifts z > 6, includ¬ 
ing the current record holder, a bright quasar at z = 7.085 
(Mortlock et al. 2011). The presence of these objects indi¬ 
cates that supermassive black holes (SMBHs) with masses of 
order 10 !) were already present in the high redshift Uni¬ 
verse at a time when the age of the Universe was only ~ 800 
million years. However, the existence of these objects is diffi¬ 
cult to understand within the standard ACDM cosmological 
model if one assumes that they grew from initial seeds with 
masses comparable to that of the Sun. Models for accretion 
onto black holes show that the characteristic growth time for 
a black hole accreting at the Eddington limit is unlikely to 
be much smaller than ~ 50 million years, unless one assumes 
an unreasonably low value for the radiative efficiency of the 
accretion process. This means that a seed black hole formed 
at very high redshift can grow by at most a factor of ~ 10' 
by the time that the high redshift SMBHs are observed. To 
produce SMBHs with the observed masses, we therefore re¬ 
quire seed black holes with masses M ~ 100 M@ or higher 
(Tanaka & Haiman 20091. In addition, if one accounts for 


the fact that radiative and mechanical feedback from the 
progenitor stars of these seed black holes will remove much 
of the gas from their local environment, it becomes difficult 
to see how the required high accretion rates can be sustained 
at high redshift, further increasing the required seed mass 
(see e.g. Alvarez, Wise & Abel 2009). 


This problem can be avoided if we assume that instead 
of forming as the remnants of Pop. Ill stars, the seed black 
holes that later grow into SMBHs form directly from the 
monolithic collapse of primordial gas. For this so-called di¬ 
rect collapse model to work, however, it is necessary that the 
collapsing gas remain warm throughout the collapse, with a 
temperature of 5000-10000 K ( Bromm fe Loeb]|2003 l. This 
keeps the Jeans mass high, preventing the gas from frag¬ 
menting into stellar mass clumps, and leading to the forma¬ 
tion of a central seed black hole with a mass of order 10 4 Mg 


or more (Begelman, Volonteri, & Rccs 2006 Choi, Shlos- 


man & Begelman 2013 Latif et al. 2013a I. For this model 


to work, we therefore require a minihalo with a virial tem¬ 
perature in excess of 10 4 K (so that the gas does not simply 
become thermally supported at low densities) that has not 
been enriched with heavy elements, which would otherwise 
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provide efficient cooling down to temperatures T <g 10 4 K 
I Omukai, Schneider & H aiman|2008 1. In addition, it is nec¬ 
essary that the formation of H 2 in the gas be strongly sup¬ 
pressed, as otherwise H 2 cooling alone would be sufficient to 


cool the gas down to temperatures far below 10 K (Omukai 


2001 Bromm & Loeb 20031. The most viable mechanism 


known for producing the required suppression of H 2 forma¬ 
tion invokes the presence of a very strong extragalactic UV 
field that immediately photodissociates almost all of the H 2 
molecules forming in the gas ( Bromm fe Loeb|2003 Visbal, 
Haiman & Bryan 20141, and that also suppresses the for¬ 
mation o f H 2 by photodissociating the in termediate H _ and 
H^ ions (Shang, Bryan fc Haiman||2010 1. 


Numerical simulations (see e.g. Bromm & Loeb 


2003 


Regan & Haehnelt 

2009; Shang, Bryan & Haiman 2010 

Latif etah] 2013a|b; 

Becerra et al. 20151 Regan, Johansson, 

& Haehnelt|2014| as well as the recent reviews by Volonteri 

2010 Haiman 2013 

and Greif 2015) have confirmed many 


of the details of this simple model, but have yet to reach 
agreement regarding the strength of the external radiation 
field that is required. It is clear that the required radiation 
field strength is orders of magnitude larger than the mean 


value produced by early generations of stars (Haiman, Abel 


& Rees 2000 Greif & Bromm 20061, but it is also clear 


that the distribution of UV radiation in the early Universe 
is highly inhomogeneous, owing to the strong clustering of 
the high-redshift protogalaxies responsible for producing it 


(Dijkstra et al. 

2008 

Ahn et al.||2009 Agarwal et al.||2012 

Dijkstra et al. 

2014 

1. The question of how rarely the re- 


quired conditions are realized in the early LTniverse there¬ 
fore depends sensitively on exactly how strong a UV field is 
actually required. This is usually quantified in terms of the 
specific intensity of the radiation field at the Lyman limit. 


Following Haiman, Abel & Rees (2000), we can measure this 

and wr ;t e it as J 21 . 


ergs ‘cm ~\ Iz ^r 


in units of 10” 

The minimum value of J 21 for which H 2 cooling is sufficiently 
suppressed is then commonly denoted as J cr it- 

Values quoted in the lite rature for J cr it range from as 
low as 20 (Inayoshi & Omukai 20111 to as high as 10 5 


(Omukai 2001). Typically, these values are determined by 


modelling the cooling and collapse of gas within minihalos 
with T v i r > 10 4 K in the presence of a strong UV background 
using 3D numerical simulations that treat the coupled ther¬ 
mal, chemical and dynamical evolution of the gas. By run¬ 
ning the simulation for a range of different values of J 21 , it is 
possible to determine the critical value at which H 2 cooling 
becomes sufficiently suppressed for the direct collapse model 
to operate. 

The scatter in the values of Jcrit obtained from these 
studies has several different causes. J cr it depends strongly 
on the spectral shape adopted for the background radia¬ 


tion field (Bromm & Loeb 2003 Shang, Bryan & Haiman 


2010) and the manner in which H 2 self-shielding is accounted 


for (Wolcott-Green & Haiman 2011), and also appears to 


vary somewhat from minihalo to minihalo (Shang, Bryan & 
Haiman|2010 Latif et al.|2014 1. However, an inportant ad¬ 
ditional contribution to the uncertainty in J cr i t comes from 
the treatment of the gas chemistry within the different stud¬ 
ies. 

In most cases, the chemical networks used in these stud¬ 
ies were originally designed for the study of Pop. Ill star 
formation in the absence of a strong extragalactic radiation 


field. Although they all include the same basic processes 
(H 2 formation via the intermediate H“ and Hi|~ ions, H 2 
destruction via collisional dissociation, charge transfer with 
H + and photodissociation, and the photodissociation of H~ 
and H^), the exact set of chemical reactions included varies 
from network to network. In addition, it is unclear whether 
any of the networks in common usage includes the full set 
of chemical processes that are important for modelling the 
formation and destruction of H 2 in an environment much 
harsher than the one that they were designed to model. 

In this paper, we attempt to identify the key reactions 
that must be accounted for when modelling the chemical 
evolution of the gas in this environment. To do this, we 
first model the chemical evolution of the gas using an exten¬ 
sive model of primordial gas chemistry that tracks 30 dif¬ 
ferent chemical species, linked by almost 400 reactions. We 
then use a reaction-based reduction technique developed by 


Wiebe, Semenov & Henning (2003) to analyze the chemistry 


and to identify the main reactions responsible for regulat¬ 
ing the H 2 abundance. This allows us to construct a much 
smaller “reduced” network that contains all of the chemical 
reactions that must be included in our chemical network if 
we are to be able to determine J cr it accurately. 

The plan of the paper is as follows. In Section [2j we 
present the details of the reaction-based reduction technique 
and the one-zone chemistry and cooling model used to sim¬ 
ulate the thermal and chemical evolution of the gas. In Sec¬ 
tion [3] we present the results of our analysis and compare 
our reduced chemical network with other chemical networks 
used in the literature to simulate the formation of black holes 
by direct collapse. Finally, we conclude in Section [4] with a 
brief summary of our main results. 


2 NUMERICAL METHOD 

2.1 Selecting the set of important reactions 

In the direct collapse scenario, the crucial quantity that de¬ 
termines whether or not the gas is able to cool significantly 
as it undergoes gravitational collapse is the abundance of 
molecular hydrogen. If the amount of H 2 that forms does not 
provide sufficient cooling, then the gas remains warm during 
the collapse, with the temperature T decreasing only slowly 


as the density increases (see e.g. Omukai 2001 Schleicher, 
Spaans & Glover 2010 k On the other hand, if enough H 2 


forms to allow the gas to cool in less than a free-fall time, 
then the temperature drops dramatically once the density 
increases above n ~ 10 2 -10 4 cm -4 . To determine the criti¬ 
cal Lyman-Werner flux, J cr it, one therefore simply needs to 
find the point at which the strength of the radiation field 
becomes large enough to suppress the H 2 abundance suffi¬ 
ciently to prevent it from cooling the gas. 

It follows from this that in order to determine J cr it ac¬ 
curately, we need to have an adequate chemical model of the 
formation and destruction of H 2 in the gravitationally col¬ 
lapsing gas. At the same time, however, we do not want to 
make our chemical model larger than is absolutely necessary, 
as each additional chemical species or chemical reaction that 
we include will further slow our numerical simulations. We 
therefore want to identify only those reactions and species 
that are most important for modelling the evolution of the 
H 2 abundance in the collapsing gas. 
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To do this, we make use of the reaction-based reduc¬ 
tion technique developed by |Wiebe, Semenov fe Henningj 
12003j) and used by them to compute the abundance of car¬ 
bon monoxide and the fractional ionization of the gas in 
simulated molecular clouds. The basic idea underlying this 
technique is quite simple. Starting from a prescribed set of 
initial conditions, we first compute the chemical and ther¬ 
mal evolution of the gas for a given value of J 21 using the 
one-zone model described in Section [2.21 below. The chem¬ 
ical network used in this one-zone model is as extensive as 
we can reasonably make it, and should therefore include all 
of the reactions that are likely to play a significant role in 
the H 2 chemistry. 

We next select our starting set of whatjWiebe, Semenov| 


& Henning (2003) term “important” species, i.e. those whose 


abundances we are interested in following accurately. In the 
present case, this set consists of only a single member - the 
H 2 molecule but in principle the same technique can be 
used with a larger set of important species. We then pro¬ 
ceed by determining the dominant production and destruc¬ 
tion reactions for each of these species at a set of different 
output times during the evolution of the gas. If the set of 
dominant reactions involve chemical species that are not in 
our starting set, then we add them to the set of important 
species, and repeat the analysis, proceeding in this fashion 
until there are no more species or reactions that need to be 
added. The result is a list of the reactions that are necessary 
at each output time in order to accurately model the abun¬ 
dances of our starting set of species. The make-up of this 
list will generally change with time: some reactions that are 
important at early times will be unimportant at late times, 
and vice versa. The final set of important reactions can then 
be obtained simply by combining the individual lists. 

To determine the dominant production and destruction 
reactions at any given output time, we proceed iteratively, 
as follows. On the first iteration, we set the weights w a of the 
important species to 1, and the weights of all other species 
to 0. We also set the weight of each reaction to zero. We then 
loop over all of the chemical species included in our chemical 
model. For each species k , we compute a new weight for each 
reaction involving that species: 


«C w (j) = max \ wr T U), ^ wT rr (k) l . (1) 

{ ^l=l,N r (k) Ul J 

Here, w a ull (k) is the current weight of species k, Wr UII (j) is 
the current weight of reaction j, iu” ew (j) is the new weight of 
reaction j, Rj is the rate per unit volume at which reaction 
j proceeds, and N r (k) is the number of reactions in which 
species k participates as either a reactant or a product. Once 
we have calculated a new weight for every reaction involving 
species k , we update its weight: 

wT"{k) = max[w“ urr (A:),max{w” ew (l),..., w*™ {N r {k))}] .(2) 


In other words, we take the new species weight to be the 
largest value from amongst the weights of the reactions in¬ 
volving that species, unless this is smaller than the current 
weight of the species. Having determined new weights for 
every species, we begin a new iteration, using the updated 
weights in Equation [l] above. We continue to iterate until 
all of the species weights have converged to within some 
pre-specified tolerance. 


One situation in which this procedure can give mislead¬ 
ing results occurs when the formation and destruction of one 
of the species is dominated by a single forwards reaction and 
its inverse. Consider, for example, the case of the formation 
of HD from H 2 

H 2 +D+ -> HD + H+. (3) 

In warm gas, this reaction occurs very rapidly, but so does 
its inverse 


HD + H^ 


H 2 + D i 


(4) 


In the extreme case in which these reactions are perfectly 
balanced, both can have large R values, and hence poten¬ 
tially also large weights, and yet the net effect of the pair 
of reactions is to leave the abundances of all of the species 
involved unchanged. In our analysis, we attempt to mitigate 
the impact of this problem by identifying reaction pairs of 
this kind and replacing the reaction rates used for these reac¬ 
tions in Equation [I] with the absolute value of the difference 
between the forward and reverse reaction rates; in the ex¬ 
ample above, this corresponds to the net rate at which H 2 
is formed or destroyed by the pair of reactions. 

Our analysis procedure leaves us with a list of weights 
for every reaction that is valid for the output time consid¬ 
ered. We can convert this into a list of reactions by retain¬ 
ing only those reactions whose weigh ts ex ceed some cut-off 
value e. As we show later in Section 


3.3 


setting e = 10 


enables us to generate a subset of reactions that allow us 
to determine J cr i t to within an accuracy of around 1% when 
compared to the results of models run with the full chemical 
network. 

By repeating the same procedure for many different out¬ 
put times, and combining the individual reaction lists, we 
can construct a reduced reaction network that is sufficient 
for modelling the chemical evolution of H 2 over the entire 
period simulated in our one-zone model. In principle, the 
reduced reaction network that we derive in this fashion is 
specific to our choice of J 21 and to the initial conditions for 
our simulation. However, by re-running the one-zone model 
and repeating the same analysis procedure for many differ¬ 
ent values of J 21 with various different sets of initial condi¬ 
tions, as described in Section [2.2| below, and the combining 
the resulting reduced networks, we can arrive at a final re¬ 
duced network that is valid over the whole range of physical 
conditions that are of interest to us. 


2.2 The one-zone model 

2.2.1 Basic details 


In order to model the chemical and thermal evolution of 
gravitationally-collapsing primordial gas illuminated by a 
strong extragalactic radiation field, we use a simple one- 
zone model that is derived from the one presented in |Glover| 
|fe Savin (|2009|. The gas density is assumed to evolve as 


d P = P_ 
d t ta ’ 


(5) 


where ts = (3-7r/32 Gp) 1 ' 2 is the free-fall time of the gas. To 
follow the thermal evolution of the gas, we write the internal 
energy equation as 
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^ = Z^ + r -A 
At p 2 At 


( 6 ) 


where e is the internal energy density, p = (7 — l)e is the 
gas pressure, V is the radiative and chemical heating rate 
per unit volume and A is the radiative and chemical cooling 
rate per unit volume. 

We compute F and A using a detailed atomic and molec¬ 
ular cooling function based on the one presented in Section 4 


of Glover & Savin (2009). In an effort to use the most up-to- 


date atomic and molecular data available, we have updated 
two of the cooling processes included in the model: the col- 
lisional excitation of H 2 by collisions with protons and by 
collisions with electrons. Full details of these updates are 
given in Appendix [A] 

To model the chemical evolution of the gas, we use a 
modified version of the chemical network present in Glover 
& Savin (2009). The original version of this network tracked 


the abundance of 30 different primordial chemical species 
and included a total of 392 reactions. We have extended this 
model by including three new reactions not treated in |Glover| 
& Savin (2009): the collisional ionization of atomic hydrogen 


by collisions with H or He, and the collisional dissociation 
of Hj by H. Details of these new reactions are given in 
Appendix [ 5 ] We have also updated several of the reaction 
rates in light of new theoretical or experimental data. Again, 
full details of these modifications are given in Appendix |B| 


2.2.2 Photochemistry and self-shielding 


The rate coefficients listed in |Glover &; S avin (2 009| ) for the 
various photochemical reactions assumed that the external 
background radiation field had the spectral shape of a 10 5 K 
black-body, cut-off at energies of 13.6 eV and above to ac¬ 
count for the effects of absorption by atomic hydrogen in 
the intergalactic medium. This choice of spectrum (which 
we refer to hereafter as a T5 spectrum, for brevity) was mo¬ 
tivated by the fact that we expect the brightest Pop. Ill stars 
to have effective temperatures of around this value (see e.g. 
|Cojazzi~et al. 2000). However, this is a reasonable approxi¬ 
mation only when the dominant contribution to the extra- 
galactic background radiation field does indeed come from 
very massive Pop. Ill stars, specffically those with masses 


greater than around 100 Mq (Bromm, Kudritzki, & Loeb 


2001). If a significant fraction of the radiation comes instead 


from less massive Pop. Ill stars or from Pop. II stars with 
systematically lower effective temperatures, then the ratio of 
the UV to the optical flux can potentially be much smaller 
than with a T5 spectrum, with important implications for 
the relative importance of H 2 photodissociation and H _ pho¬ 
todetachment. Several authors have attempted to quantify 
the dependence of J cr it on the spectral shape of the extra- 
galactic radiation field by performing simulations both with 
a T5 spectrum and with a T4 spectrum: a 10 4 K diluted 
black-body with a similar high-energy cut-off to the T5 spec¬ 
trum. These studies (see e.g. Bromm & Loeb 2003; Shang, 


Bryan fc Haiman|2010 1 find that the choice of spectrum has 


a large influence on the value of J cr it ■ For instance, |Shang, 


Bryan & Haiman (2010) find that in their models, the value 
of Jcrit for a T5 spectrum lies in the range 10 4 < J cr it < 10 s , 
but that for a T4 spectrum, 30 < J cr it < 300, a difference of 
two orders of magnitude or more. 


To account for this uncertainty, we perform two sets of 
simulations: one set using a T5 spectrum and a second set 
using a T4 spectrum. For the runs with the T5 spectrum, 
we use the photochemical rates listed in |Glover fe Savin] 
12009) for all processes except for the photodissociation of 
I l(T. which is treated with the density-dependent approach 
described in Appendix |B2| For the other set of runs, we have 
recomputed the photochemical rates using a T4 spectrum 


and the same basic approach as in Glover & Savin] (2009). 


In practice, of course, both the T4 and T5 spectra 
are relatively crude approximations to the actual spectrum 
of the extragalactic background produced by high-redshift 
star-forming protogalaxies (|Agarwal fc~Khochfar 2015] [Sug7| 


imura, Omukai & Inoue 2014 Latif et al. 2015). However, 


they should bracket the behaviour seen in more realistic 
models, making them suitable choices for the purposes of 
our study. 

When treating most of the photochemical reactions, we 
assume that the gas remains optically thin, as the continuum 
opacity of low density primordial gas is very small (see e.g. 
Lenzuni, Chernoff & Sal peter||199 1). The important excep¬ 
tion is H 2 photodissociation, as the rate of this processes can 
be strongly affected by H 2 self-shielding. To account for this, 
we use the modified form of the 


Draine & Bertoldi 


(1996) 


shielding function introduced by | Wolcott-Green, Haiman &| 
Bryan] |2011): 


f s h(Nn 2 ,T ) = 


0.965 


+ ■ 


0.035 


( l+x/b 5 ) 14 (1 + x ) 0 - 5 

x exp [—8.5 x 10 _4 (1 + x) 0 ' 5 ] , 


( 7 ) 


where x = Nh 2 /5 x 10 14 cm~ 2 , 65 = 6/10 5 ems -1 , and b 
is the Doppler broadening parameter, which we assume to 
be dominated by the effects of thermal broadening. In or¬ 
der to evaluate this expression for / a h, we need to know the 
column density of H 2 in the collapsing gas cloud. We com¬ 


pute this using a similar method to Omukai (2001): we as¬ 


sume that the dominant contribution to the shielding comes 
from a core region with radius R c = A.j, where Aj is the 
Jeans length. The required H 2 column density then follows 
as A r H 2 = nn 2 R c = 71 h 2 Aj, where uh 2 is the number density 
of H 2 in our one-zone model. 

It should be noted that the simple approximation that 
we use here to compute Nh 2 has been shown to overestimate 
the actual amount of shielding by up to an order of magni¬ 
tude in some cases compared to the results of fully 3D treat¬ 
ments (Wolcott-Green, Haiman & Br yan|2 011). For this rea¬ 
son, the absolute values of J cr i t that we derive in this study 
should be treated with considerable caution. However, this 
simplification should not affect our conclusions regarding 
the set of chemical reactions that are required to accurately 
determine J cr it, or our findings regarding the sensitivity of 
Jcrit to uncertainties in the rates of these reactions. 


2.2.3 Initial conditions 

We perform simulations using three different set of initial 
conditions, as detailed in Table [l] In runs 1 and 4, we start 
with a low initial temperature, To = 200 K, and an initial 
fractional ionization and H 2 fractional abundance close to 
those found in the IGM prior to the onset of Pop. Ill star 
formation (a; e ,o = 2 x 10 4 , :eh 2 ,o = 2 x 10~ 6 ). In runs 2 and 
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Table 1. List of simulations 


Run 

To (K) 

x e,0 

*^H2 ,0 

Spectrum 

1 

200 

2 x 10- 4 

2 x 10 -6 

T4 

2 

8000 

2 x 10“ 4 

2 x 10 -6 

T4 

3 

8000 

1.0 

0.0 

T4 

4 

200 

2 x 10 -4 

2 x 10 -6 

T5 

5 

8000 

2 x 10“ 4 

2 x 10- 6 

T5 

6 

8000 

1.0 

0.0 

T5 


5, we start with the same chemical composition but with a 
much higher gas temperature, To = 8000 K. Finally, in runs 
3 and 6, we set To = 8000 K, but assume that the gas was 
initially fully ionized (x e ,o = 1.0, xh 2 ,o = 0.0). The initial 
density in all of the simulations is set to no = 0.3 cm -3 , cor¬ 
responding approximately to the mean density within the 
virial radius of a protogalaxy forming at a redshift z = 20. 
We have explored the effects of adopting a larger initial den¬ 
sity but find that this makes little difference to our results. 

In all of our simulations, we adopt an elemental abun¬ 
dance (relative to hydrogen) of An = 2.6 x 10 -5 for deu- 
terium and Am = 4.3 x 10 -10 for lithium ( Cyburt|2004 1. We 
set the initial abundances of D + and HD so that they are a 
factor of An smaller than the initial H + and H 2 abundances, 
respectively. We assume that the lithium was initially en¬ 
tirely in neutral atomic form, and set the initial abundances 
of all of the other chemical species to zero. 

For each set of initial conditions, we run two different 
sets of simulations, one with a T4 spectrum and the other 
with a T5 spectrum, as indicated in Table [I] 


3 CONSTRUCTION OF AN ACCURATE 
REDUCED NETWORK 


3.1 Determination of J cr it 


Before we can identify the set of reactions that it is neces¬ 
sary to include in the chemical model in order to accurately 
determine J cr it, it is first necessary to establish the value of 
J crit for each of our six different setups. We do this using a 
binary search method. We start by selecting two values of 
J 21 that are certain to bound J cr it: a low value, JA.iow = 1, 
that we have verified is too small to prevent efficient H 2 cool¬ 
ing from occurring in any of the models, and a high value, 
Jn.high = 10 4 , that completely suppresses H 2 cooling in all 
of the models. We then compute a new value of J 21 using 
the equation 


J21,new — (J21,low X Til,high) , 


( 8 ) 


and run a simulation using this new value. If Tgnew is large 
enough to prevent efficient H 2 cooling, then we adopt it as 
our new value of Ti.high; otherwise, we take it as our new 
value of J 2 i,iow- We proceed in this fashion until the dif¬ 
ference between Tpiow and Ti.high is less than 0.2%, and 
take the final value of Jn.uew as our estimate for Jcrit■ The 
resulting values are listed in Table [2] for each of our six dif¬ 
ferent setups. In agreement with previous work, we find that 
the choice of spectral shape has a large influence on J cr it- 
With a T4 spectrum, we find that J cr it ~ 17-18, while a T5 
spectrum yields J cr it ~ 1630. We see also that the results 


Table 2. Critical Lyman-Werner flux as a function of e 


Run 

£ = 0 

1 

O 

II 

«/crit 

e = 10 -3 

e = 0.01 

e = 0.1 

1 

17.0 

17.1 

17.1 

17.1 

17.5 

2 

18.0 

18.0 

18.0 

18.1 

18.6 

3 

18.0 

18.1 

18.1 

18.1 

18.6 

4 

1640 

1640 

1640 

1650 

4260 

5 

1630 

1630 

1630 

1640 

4260 

6 

1630 

1630 

1630 

1640 

4260 


The values shown here are specified to only three significant 
figures, and were computed using chemical networks consisting 
of all chemical reactions with maximum weights exceeding the 
listed value of e. The case e = 0 corresponds to the full chemical 
model. 


that we obtain are not particularly sensitive to our choice of 
initial conditions. 

It is interesting to compare our results to those of previ¬ 
ous one-zone studies. For the T4 spectrum, our finding that 
Jcrit ~ 17-18 is in fairly good agreement with the values 


Jcrit = 20,25,39 derived by Inayoshi & Omukai (20111, Sug 


imura, Omukai & Inoue (20141, and Shang, Bryan & Haiman 


(2010), respectively. The difference between our value and 


the values derived in these other studies is most likely due 
to the impact of the differences in the chemical model and 
the choices made for some of the key rate coefficients, since 


as we will see later in this paper and in Paper II (Glover 


20151, these differences can easily introduce an uncertainty 


of a factor of a few into our estimates for J cr it- 

For the T5 spectrum, our value of J cr it is almost an or¬ 
der of magnitude smaller than th e values of Jcrit = 1.2 x 10 4 
and 1.6 x 10 4 found by Shang, Bryan & Haiman (2010) and 


Inayoshi & Omukai (20111, respectively. However, this dif¬ 


ference can be explained by differences in our treatment of 
H 2 self-shielding. Both of these previous studies used the 
H 2 self-shielding function computed by |Draine fc Bertoldi| 
(|1996[), while we used instead the modified version given in 


Wolcott-Green, Haiman & Bryan (2011), which more accu¬ 


rately treats the self-shielding of H 2 in hot gas. |Sugimurig1 
|Omukai &; Inou c (2014) performed one-zone calculations 
with a T5 spectrum using both of these treatments of H 2 
self-shielding, and showed that the value of J cr ;t that they 
obtained with the Wolcott-Green, Haiman & Br yan| pOll) 
treatment was approximately an order of magnitude smaller 
than that obtained using the Draine & Bertoldi ( 1996|) treat- 
ment. The value they obtained for J cr it with the |Wolcott-| 


Green, Haiman & Bryan (2011) treatment was J cr it — 1400, 


in good agreement with the value we find using our model. 


3.2 The reduced network 

Having determined the value of Jcrit in each run, we next 
analyze the full set of chemical reactions taking place during 
the evolution of the gas, as outlined in Section [2.1| For each 
run, we consider a large set of output times, and use our 
computed values for the density, temperature and chemical 
composition of the gas to determine the weight of each re¬ 
action in our full set at that output time. We record these 
weights, and repeat this procedure for a large number of 
different output times. We then determine the maximum 
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weight for each reaction in our model. To help ensure that 
our reduced network will be robust against minor changes in 
the physical conditions in the gas, we consider not only the 
case when J 21 = Jcrit, but also perform the same analysis 
for simulations with J 21 = 0.3J cr it and J 21 = 3J cr it, tak¬ 
ing the maximum weight for each reaction to be the largest 
of the weights obtained for the reaction with these three 
different setups. Finally, we construct our reduced network 
using only those reactions whose maximum weights exceed 
e = 10 -4 . The resulting set of reactions for each run is listed 
in Table [3] along with the corresponding set of maximum 
weights. 

We see immediately that there is a subset of reactions 
with very large weights that appear in the reduced network 
for each of our runs. This is not surprising: most of the reac¬ 
tions in this set play such important roles in the regulation 
of the H 2 abundance that any model omitting them could 
not hope to be representative of the real behaviour of the 
gas. Examples include the photodissociation and collisional 
dissociation of H 2 , the formation of H 2 from H“, or the for¬ 
mation and photodissociation of the H _ ion itself. 

There are also a few reactions that have relatively high 
weights in some of our runs that are more unexpected. Per¬ 
haps the best example of these is the collisional ionization 
of atomic hydrogen by collisions with other hydrogen atoms, 
i.e. 

H + H -> H + H+ +e“. (9) 

Most previous networks developed to treat primordial chem¬ 
istry have omitted this reactiorQ assuming it to be unimpor¬ 
tant in comparison to the collisional ionization of hydrogen 
by electrons, since the latter reaction has a far larger rate 
coefficient. However, while it is certainly true that collisions 
with electrons dominate in many circumstances - for exam¬ 
ple in gas cooling and recombining from an initially ionized 
state - this is not true in all of the runs we examine here. 
In particular, if the initial electron abundance is small, as 
for example in runs 1 or 2, then collisions between hydrogen 
atoms happen so much more frequently than collisions be¬ 
tween hydrogen atoms and electrons that this reaction can 
come to dominate the ionization rate despite its small rate 
coefficient. 

To illustrate the importance of this reaction, we re-ran 
our determination of J cr i t for models 1 and 4 with the value 
of its rate coefficient set to zero. We found that in this case 
Jcrit = 11.8 for run 1 and J cr i t = 1280 for run 4. In other 
words, omitting this reaction leads to a 20-30% decrease in 
the derived value of J cr it when we start from low ionization 
initial conditions. If we also disable the reaction 

H + He ->• H + + e" + He, (10) 

which again is often omitted from treatments of primordial 
chemistry, then J cr it is reduced even further, to J cr it = 8.9 
or 1110 for runs 1 and 4, respectively. 

We can understand why these two reactions appear to 
be so important for determining J cr it if we examine how 

1 The notable exception is jOmuka i] (200 1 }. who does include this 
process in his chemical network. Studies of SMBH formation that 
make use of his network (e.g. |Sugimura, Omukai A I noue| |20f4| 
therefore already account for this reaction. 



Figure 1. Top panel: evolution of the gas temperature as a func¬ 
tion of the hydrogen nuclei number density in runs with the same 
setup as in run 6 (a T5 spectrum, fully ionized initial conditions). 
Results are shown for J 21 = 1630 (solid line) and J 21 = 1635 
(dashed line). Bottom panel: evolution of the fractional ioniza¬ 
tion as a function of density in the same runs. 

the fractional ionization of the gas evolves as a function of 
density when J21 is close to Jcrit • In Figure [l] we show how 
the temperature and the fractional ionization evolve with 
density in two simulations carried out with the same initial 
conditions as run 6 (i.e. fully ionized gas, illuminated by a 
T5 spectrum), with J21 = 1630 (dashed line) and J21 = 
1635 (solid line). These two values of J21 were chosen to 
tightly bracket J cr it- We see that in both simulations, the 
temperature evolution is essentially identical until we reach 
a density of n ~ 10 3 cm~ 3 . At this point, the behaviour of 
the two runs diverges. In the run with J21 < Jcrit, the H2 
fraction in the gas at this point is high enough to allow it 
to begin to cool significantly. As it cools, the effects of H 2 
collisional dissociation become less effective, allowing the H 2 
fraction to increase further, and so the cooling runs away, 
rapidly driving the temperature down to T ~ 1000 K. In 
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Table 3. List of reactions with maximum reaction weights greater than 10 4 in at least one simulation 


No. 

Reaction 

Run 1 

Run 2 

Maximum reaction weight 
Run 3 Run 4 

Run 5 

Run 6 

1 

H 2 + 7 -4- H + H 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

2 

H 2 + H-S.H + H + H 

0.49 

0.50 

0.49 

0.49 

0.49 

0.49 

3 

H” + H —> H 2 + e” 

0.47 

0.47 

0.47 

0.50 

0.47 

0.47 

4 

H+ + H -4 H 2 + H+ 

0.41 

0.49 

0.48 

4.7 x 10” 2 

4.3 x 10” 2 

4.4 x 10” 2 

5 

H + + e — —> H + 7 

0.27 

0.12 

0.48 

4.7 x 10” 3 

1.0 x 10” 2 

4.3 x 10” 2 

6 

H + e — —»■ H _ + 7 

0.23 

0.23 

0.23 

0.25 

0.23 

0.24 

7 

H — + 7 —»■ H + e — 

0.21 

0.15 

0.14 

0.25 

0.23 

0.23 

8 

H + H+ -x H+ + 7 

0.20 

0.24 

0.24 

2.3 x 10” 2 

2.2 x 10” 2 

2.1 x 10” 2 

9 

H+ + 7 -> H+ + H 

0.12 

0.24 

0.24 

7.2 x 10” 3 

2.1 x 10” 2 

2.0 x 10” 2 

10 

H + H—>H++e”+H 

9.5 x 10” 2 

0.12 

1.1 x 10” 2 

1.4 x 10” 2 

1.1 x 10” 2 

1.4 x 10” 3 

11 

H-+H-)H + H + e" 

8.8 x 10” 2 

8.9 x 10” 2 

8.8 x 10” 2 

9.3 x 10” 2 

9.2 x 10” 2 

9.2 x 10” 2 

12 

H + e” -7 H+ + e” + e” 

6.1 x 10” 2 

0.22 

4.0 x 10” 2 

7.8 x 10” 3 

2.0 x 10” 2 

3.9 x 10” 3 

13 

H+ + He -» HeH+ + H 

3.5 x 10” 2 

2.8 x 10” 3 

2.7 x 10” 3 

3.6 x 10” 4 

3.0 x 10” 4 

3.0 x 10” 4 

14 

H + He — ► H+ + e” + He 

2.7 x 10” 2 

3.6 x 10” 2 

3.1 x 10” 3 

3.3 x 10” 3 

3.2 x 10” 3 

3.9 x 10” 4 

15 

H 2 + H+ — > H+ + H 

1.0 x 10” 2 

8.0 x 10” 3 

4.2 x 10” 2 

1.1 x 10” 3 

1.1 x 10” 3 

1.1 x 10” 3 

16 

H 2 + He -> H + H + He 

5.9 x 10” 3 

5.9 x 10” 3 

5.9 x 10” 3 

5.8 x 10” 3 

5.8 x 10” 3 

5.8 x 10” 3 

17 

HeH+ + H — t H+ + He 

2.3 x 10” 3 

2.8 x 10” 3 

2.6 x 10” 3 

3.6 x 10” 4 

3.0 x 10” 4 

3.0 x 10” 4 

18 

H + H + H-xHa + H 

1.4 x 10” 3 

1.3 x 10” 3 

1.3 x 10” 3 

— 

— 

— 

19 

H” + He -> H + He + e” 

1.3 x 10” 3 

1.4 x 10” 3 

1.3 x 10” 3 

2.5 x 10” 3 

2.0 x 10” 3 

1.9 x 10” 3 

20 

H++H-S-H + H++H 

1.3 x 10” 3 

3.6 x 10” 4 

5.2 x 10” 4 

— 

— 

— 

21 

He + H+ — ► HeH+ + 7 

5.3 x 10” 4 

— 

— 

— 

— 

— 

22 

H” + H+ -4 H + H 

4.6 x 10” 4 

4.6 x 10” 4 

4.0 x 10” 4 

7.5 x 10” 4 

1.4 x 10” 3 

6.6 x 10” 2 

23 

H+ + e” —t H + H 

1.7 x 10” 4 

2.0 x 10” 4 

2.4 x 10” 2 

— 

— 

8.7 x 10” 4 

24 

HeH+ + e” ->• He + H 

— 

— 

1.9 x 10” 4 

— 

— 

— 

25 

H” +H+ -4 H+ +e” 

— 

— 

1.0 x 10” 4 

— 

— 

9.6 x 10” 4 

26 

H” + e” —X H + e” + e” 

— 

— 

— 

1.0 x 10” 4 

2.6 x 10” 4 

9.5 x 10” 3 


The reactions with numbers listed in bold font make up the minimal reduced model described in Section |3.3| 


the run with J 21 > Jcrit, on the other hand, the smaller H 2 
abundance means that H 2 cooling never becomes effective 
and the gas remains hot throughout the simulation. 

The value of J cr i t is therefore determined in this case 
by the chemical state of the gas at n ~ 10 3 cm -3 and 
T ~ 7500 K. The H 2 fraction in the gas in these condi¬ 
tions depends on the fractional ionization. Figure [l] shows 
that even though the gas is initially fully ionized, by the 
time we reach n ~ 10 3 cm -3 , the fractional ionization has 
lost its memory of this initial state and has decreased to 
i~4x 10” 5 . This value is set by the balance between the 
collisional ionization of hydrogen by collisions with electrons, 
H atoms and He atoms (i.e. reactions 12, 10, and 14, respec¬ 
tively) and the radiative recombination of H + . If we compare 
the rates of reactions 10, 12 and 14 in these conditions, we 
find that Riq ~ 4x 10” 16 cm” 3 s” 1 , R 12 — 2x 10” 16 cm” 3 s” 1 
and i?i 4 ~ 10” lb cm” 3 s” 1 . In other words, because of the 
very low fractional ionization of the gas, reactions 10 and 
14 between them provide around 70% of the total collisional 
ionization rate, with the main contribution coming from re¬ 
action 10. Omitting these reactions therefore leads to a fac¬ 
tor of ~ 3 error in the total ionization rate, and hence a 
factor of ~ 2 error in the fractional ionization. This reduces 
the H 2 formation rate by a similar amount and hence also 
reduces our estimate of J cr it- 

Another interesting thing that Table [3] shows us is that 
once we start to look at reactions with lower weights, which 
are less important overall for our determination of J cr it, we 
see substantially more variation from run to run. Only for 
runs 4 and 5 does our reaction weighting algorithm provide 


us with exactly the same set of reactions with maximum 
weights greater than e; for the other runs, we see minor 
differences, generally involving a small number of processes 
with weights close to our cut-off. In order to have a reaction 
network that is robust to changes in the initial conditions 
and details of the background radiation field, we therefore 
recommend that one uses all of the 26 reactions listed in 
Table H 


3.3 Effect of varying e 

In order to establish that our choice of a value of e = 10” 4 
for the reaction weight cut-off in our selection algorithm 
does indeed allow us to select all of the reactions important 
for determining J cr it, we re-ran our models using a reduced 
chemical network consisting only of the 26 reactions listed 
in Table [3] and recalculated J cr it for each model using the 
same approach as before. The resulting values are shown 
in Table [2] We see that there are minor differences in the 
predicted values of J cr it in a couple of the runs, but that they 
are at the level of less than 1%. This is much smaller than 
the scatter in J cr i t from halo to halo caused by differences in 
the dynamical evolution of the gas ( Latif et al.| [2014). This 
demonstrates that our 26 reaction reduced network is indeed 
accurate enough for our purposes. 

The question then naturally arises as to whether we can 
simplify our reduced network even more. If we increase the 
value of e, and hence select even fewer reactions to join the 
set of “important” reactions, then how badly does this com¬ 
promise our ability to determine Jcrit accurately? To explore 
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this, we constructed several additional reduced networks, 
in which we retained only those reactions with maximum 
weights greater than e = 10 —3 , 0.01 or 0.1 in at least one 
run, and then used these even more highly simplified net¬ 
works to determine J cr it- When performing this analysis for 
the t = 0.01 case, we found that in order to produce mean¬ 
ingful results, it was necessary to include one reaction - the 
destruction of HeH + by collisions with H, reaction 17 in the 
Table - that formally has a weight that is less than 0.01. 
The reason for this is that otherwise our e = 0.01 model 
would contain a reaction that forms HeH + (reaction 13 in 
the Table), but no reaction that destroys it, meaning that 
the abundance of HeH + would increase indefinitely. 

The results that we obtained when varying t are shown 
in columns 4-6 of Table [2] We see that the results we obtain 
for e = 10~ 3 or £ = 0.01 are very similar to those we obtain 
in the e = IIP 4 case, with the derived values of J cr it differing 
by no more than around 1%. However, in the e = 0.1 case, 
we see a clear change in behaviour: the error in the value 
of Jcrit we obtain for the T5 spectrum has increased from 
~ 1% to a factor of a few. This suggests that if we want to 
simplify the chemistry as much as possible while still being 
able to determine J cr it accurately, then the best choice is 
a minimal reduced network consisting of reactions 1-15, 17, 
22, and 23. On the other hand, if we want to be a little more 
conservative, in view of the fact that the physical conditions 
encountered in more realistic simulations are not perfectly 
reproduced by our one-zone model, then we should retain 
all 26 of the reactions listed in Table [3] 


Table 4. Dependence of ,/ ( ,rii. on our treatment of [if excitation 


Run 

Fiducial 

Jcrit 
v = 0 

LTE 

1 

17.0 

17.6 

14.2 

2 

18.0 

18.6 

14.5 

3 

18.0 

18.6 

14.5 

4 

1640 

1650 

1560 

5 

1630 

1640 

1560 

6 

1630 

1640 

1560 


comparison. We see that there is only a small difference be¬ 
tween the values of J cr it that we obtain using our fiducial 
model and those from the v = 0 model. This is not surpris¬ 
ing, since the question of whether or not enough H 2 forms to 
cool the gas is largely determined by the chemistry occurring 
at densities n < n cr i t , where the v = 0 rates are a good ap¬ 
proximation. If we instead adopt LTE rates throughout, we 
find a larger difference between the models, approximately 
20% in the runs with the T4 spectrum and ~ 5% in the runs 
with the T5 spectrum. As this represents a worst case, the 
real error introduced is almost certainly smaller than this. 
Our simplified treatment of excitation therefore does 
not represent a major source of error in our determination 
Of Jcrit- 


3.4 Sensitivity to the treatment of excitation 


3.5 Importance of dissociative tunneling 


As we discuss in more detail in the Appendix, the rate at 
which Hi!" is destroyed by processes such as dissociative re¬ 
combination or photodissociation is highly sensitive to the 
vibrational state of the molecular ions: it is far easier to de¬ 
stroy Hj ions in states with d>0 than ones that are in the 
v = 0 vibrational ground state. Therefore, any model of the 
chemistry of Hj in primordial gas should also account, at 
least approximately, for the degree of excitation of the ions. 
In the chemical model presented in this paper, this is done in 
a simple but approximate fashion, for reasons of computa¬ 
tional convenience. However, in principle one can construct 


far more complex and accurate models (see e.g. Coppola 
et al. 2011). An obvious question is therefore whether this 


complexity is necessary if one wants to accurately determine 
Jcrit , or whether the simple treatment used in this paper suf¬ 
fices. 

To answer this question, we performed two additional 
sets of runs. In the first set, we artificially set the H^" critical 
density to a very large number, so that all of the chemical 
rates used for the ion were the ones for Hj in its ground 
state. In the second set of runs, we instead set the H % critical 
density to zero, so that the local thermodynamic equilibrium 
(LTE) rates were used for all of the reactions involving H(J" as 
a reactant. These two limiting cases should bracket the real 
behaviour of the Hj ions in the gas, and so by comparing 
the results of these two sets of runs, we can establish how 
sensitive J cr it is to the details of the excitation of the Hj 
ions. 

We show the results from these two sets of runs in Ta¬ 
ble [4] along with the results from our fiducial model for 


Our analysis in Section [3.2| showed that the collisional dis¬ 
sociation of H 2 has one of the largest weights in our reduced 
chemical network, and therefore it is very important to treat 
this process accurately. As discussed at some length in e.g. 


Martin, Schwarz & Mandy (1996), there are two distinct 


processes that contribute to the total collisional dissocia¬ 
tion rate. The first is direct collisional dissociation, where 
the H 2 molecule undergoes a transition from a bound state 
directly into the continuum of classically unbound states. 
The second is dissociative tunneling, where the H 2 molecule 
is excited into a quasi-bound state - a state that has an 
internal energy larger than is required for dissociation, but 
which is separated from the continuum by a barrier in the 
effective potential. Quantum tunneling through this barrier 
then leads to spontaneous dissociation. In low density gas 
this process of dissociative tunneling can make a major con¬ 
tribution to the total collisional dissociation rate. 

Some studies of the direct collapse model for the forma¬ 
tion of massive black holes include only the effects of direct 
collisional dissociation, and not those of dissociative tunnel¬ 
ing (e.g. Shang, Bryan & Haiman 2010), but as Latif et al. 
12014) point out, this potentially introduces significant un¬ 
certainty into the resulting values derived for Jcrit ■ In an 
effort to quantify this uncertainty, we re-ran our set of six 
simulations using our reduced chemical model, but with the 
effects of dissociative tunneling disabled. We found that in 
this case, J cr it = 33.1 in runs 1 -3 and J cr it = 2780 in runs 
4-6. In other words, the failure to account for this process 
when modelling H 2 dissociation in this scenario can intro¬ 
duce an error of almost a factor of two into J cr it. 
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Table 5. Comparison of J cr it determined using our reduced 
model with that determined using several other simplified models 
from the literature 


networks to consider: the react-primordial_photoH2 network 
and the react-xrays network. 

The react-primordiaLphotoH2 network consists of reac- 


Run 

This work 

KROME 

'^crit 

ENZO 

Omukai (2001) 

network, along with several additional reactions, largely in¬ 
volving the ionization and recombination of helium and the 

1 

17.1 

8.9 

21.5 

17.1 

chemistry of D, D + and HD, that our results demonstrate 

2 

18.0 

14.7 

26.7 

18.3 

are not important for the determination of J cr it- The re- 

3 

18.1 

14.9 

26.8 

18.3 

act-xrays network consists of the reactions mentioned above, 

4 

1640 

1120 

1930 

3040 

plus two more from our reduced network, reactions 9 and 18. 

5 

1630 

1160 

1940 

3040 

It also accounts for charge transfer between hydrogen and 

6 

1630 

1160 

1940 

3040 

helium, and includes a more extensive treatment of deu- 

3.6 

Comparison with other chemical models 

terium chemistry, but once again, our results suggest that 
these processes are unimportant for determining J cr it- 

We therefore see that compared to our reduced network, 


It is interesting to investigate whether there are any sig¬ 
nificant differences between the reduced chemical network 
constructed in this paper and the other chemical networks 
in common usage in studies of the direct collapse model 
for the formation of massive black holes. The majority of 
existing studies use one of three treatments: the krome as- 
trochemistry package ( Grassi et al.|2014 |, which provides a 
range of different networks for modelling primordial chem¬ 
istry; the primordial chemistry network implemented in the 
enzo adaptive mesh refinement code, described in |Bryan| 
et al. (20141, which is a modified and extended version of 
the network outlined i n | Abel et al.| ( |1997| ); or the network 
introduced in[Omukai (20011. We compare our reduced net¬ 
work to each of these treatments below. 


3.6.1 The krome package 

The krome astrochemistry package is a Python-based pre¬ 
processing system that converts a simple textual description 
of an astrochemical network into a set of subroutines for the 
solution of the chemical rate equations for the species con¬ 
tained in that network. In addition, krome also contains 
support for a large number of different heating and cooling 
processes, allowing the thermal energy equation to be solved 
along with the chemical rate equations, krome is still un¬ 
dergoing active development, but in our discussion here we 
refer to the state of the code at the time of writing^] note 
that this differs in some respects from the previous version 
of the package described in Grassi et al. ( 2014| . 

The krome package provides a number of example net¬ 
works, including several designed for modelling primordial 
gas. Unfortunately, these come with very little documenta¬ 
tion, and so it is not immediately obvious which network is 
the best choice for which application, or even which (if any) 
of these networks has been used in studies such as ILatifl 


et al. (20141, Van Borm et al. (20141, or Latif et al. (2015) 


that have been carried out using krome. However, some in¬ 
vestigation allows us to immediately cut down the number 
of possibilities. We can clearly ignore any of the primordial 
networks that do not include the photodissociation of H~ 
or H 2 , since these are crucial processes in the case of the 
T4 or T5 spectra, respectively. This leaves us with only two 


2 Specifically, on January 3rd, 2015, at which time the most re¬ 
cent commit was 53a32c2ab5d53994f e4893676dbfbcf 08125c23d 


the react-xrays network is missing reactions 10, 13-14, 16- 
17, 19-21 and 24, while the react_primordiaLphot.oH2 net¬ 
work is missing these plus also reactions 9 and 18. To quan¬ 
tify the effect of omitting these reactions on the determina¬ 
tion of Jcrit, we have rerun our set of six simulations using 
only the reactions contained in the react-xrays network, and 
determined the value of J cr it in each case. Note that when 
doing so, we use the rate coefficients for each reaction taken 
from our chemical model, which are not always the same 
as those adopted by the krome model. This is because we 
are interested here only in the effects of omitting some of 
the reactions that we have determined are important, and 
not on quantifying the uncertainty introduced into J cr it by 
differences in our choice of rate coefficients. We examine the 
latter issue at some length in a companion paper (jGlover 


20151. 


The results of our comparison are shown in Table [5] We 
see that the krome network systematically underestimates 
J cr it. The discrepancy is worst in run 1, where the error is 
almost a factor of two, but even in the best case, the error 
is ~ 20%. We have also carried out a similar comparison 
using the set of reactions in the react_primordial_photoH2 
network, but find that in this case we recover very similar 
results. We have investigated the source of this discrepancy 
and find that almost all of it is accounted for by the omission 
of reactions 10 and 14 from the krome model. If we add 
these two reactions to the set included in the react-xrays 
network, then we can reproduce the values of Jcrit that we 
obtain with our minimal model to within ~1%. 

3.6.2 The enzo network 

The primordial chemistry network implemented in the enzo 
code consists of the same subset of the reactions from our 
minimal model as in the krome react-xrays network, i.e. re¬ 
actions 1-8, 11, 12, 15, 22, 23, 25 and 26. As in the case of the 
krome networks, it also includes a number of reactions in¬ 
volving the ionization and recombination of helium and the 
chemistry of deuterium that play no role in the determina¬ 
tion of Jcrit- The only difference between the enzo network 
and the the krome react-xrays network that is relevant for 
our present study is that the krome network accounts for 
the collisional dissociation of H 2 via dissociative tunneling 
as well as direct dissociation into the continuum, whereas 
the enzo model only includes the latter process. 

In Table [5] we compare the values of J cr it that we ob¬ 
tain using the ENZO model with those that we obtain using 
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our reduced model. We find that the enzo model systemati¬ 
cally overestimates J cr it by between 20% and 50%. Interest¬ 
ingly, the mean difference is smaller than with the krome 
model, despite the fact that the enzo model is less com¬ 
plete than the krome model. The reason for this is that 
while the neglect of reactions 10 and 14 in the ENZO model 
tends to decrease our estimate of J cr it, the neglect of dis¬ 
sociative tunneling has the opposite effect, as we have seen 
already in Section |3.5| Therefore, the two effects cancel to 
some extent, and so the overall difference with the results of 
our reduced model is less than we would initially expect. 


3.6.3 The Omukai (2001) network 


The chemical network adopted by Omukai (2001) in his pi¬ 


oneering study of the effects of a strong extragalactic ra¬ 
diation field on the gravitational collapse of primordial gas 
includes reactions 1 10, 12, 15, 18, 22 and 26 from our re¬ 
duced network. It also includes a number of additional re¬ 
actions such as the three-body recombination of H + , or the 
collisional ionization of electronically excited hydrogen that 
are unimportant at the densities at which the value of J cr it is 
set, but which play important roles in regulating the ioniza¬ 
tion state and thermal evolution of the gas at much higher 


densities. The Omukai (2001) network, which for the sake of 


brevity we refer to hereafter as the O01 network, or an up¬ 
dated version of it, has subsequently been used in a number 
of different studies of aspects of the direct collapse model 


for black hole formation (e.g. 

Omukai, Schneider & Haiman 

2008{ Inayoshi & Omukai 2011 2012 Inyoshi, Omukai & 

Tasker|2014 Sugimura, Omukai & Inoue|2014 

)- 


when we use the set of reactions contained in the O01 net¬ 
work with the values obtained using our reduced network. 
In performing this comparison, we have assumed that the 
O01 network accounts for both the direct collisional disso¬ 
ciation of H 2 and its destruction by dissociative tunneling. 
This was not the case in the original version of the network, 
which used a highly outdated treatment of H 2 collisional 
dissociation taken from Lepp & Shull ( 1983|. More recent 


studies use an updated treatment based on |Martin, Schwarz| 
| fe Man dy ( 1996]) , but do not clarify whether they include 
only the direct collisional dissociation or also the dissocia¬ 
tive tunneling term. 

We see from Table [5] that there is very good agreement 
between the results of our reduced network and those ob¬ 
tained using the O01 network in the case of the simulations 
performed using the T4 spectrum. In these runs, the maxi¬ 
mum difference between the two models is around 2%. The 
main reason that the O01 model performs much better than 
the krome or enzo models in this case is that it includes 
the effects of reaction 10, which, as we have already seen, 
significantly contributes to the ionization state of the gas at 
the densities where the value of J cr it is set. 

I 11 the case of the runs performed using the T5 spec¬ 
trum, however, we see that there is a difference of almost a 
factor of two between the results from our reduced model 
and those from the O01 model. This difference is driven 
almost entirely by the fact that the O01 model does not 
include reaction 11 from our reduced model, the collisional 
detachment of H - by H: 


H~+H—s-H + H + e". (11) 

Although this reaction generally occurs more slowly than 
the associative detachment reaction responsible for forming 
H 2 from H - (reaction number 3 in our reduced network), 
at the temperatures reached in gas with J 21 ~ Jcrit, the dif¬ 
ference between the two rates is less than an order of mag¬ 
nitude. In the runs with the T4 spectrum, both reactions 
have rates that are small compared to the photodetachment 
rate (reaction 7), and hence including reaction 11 in the 
model makes little difference to the equilibrium H~ abun¬ 
dance or the H 2 formation rate. On the other hand, in runs 
performed using the T5 spectrum, H~ photodetachment is 
much less important and so reaction 11 plays a much more 
important role in regulating the H“ abundance. Its omission 
from the chemical network leads to one overestimating the 
H~ abundance and hence overestimating the H 2 formation 
rate. Consequently, a larger value of J21 is required in or¬ 
der to suppress H 2 cooling, and so the resulting values for 
Jcrit are systematically larger. We have verified this by re¬ 
running the T5 models with a modified version of the O01 
network that includes reaction 11. We find that in this case, 
the difference between the value of J cr it that we obtain for 
the modified O01 model and the one that we obtain for our 
reduced model is only ~5%, with this remaining difference 
mostly likely being due to the fact that the O01 model does 
not include the effects of the collisional ionization of H by 
collisions with He (reaction 14). 


4 SUMMARY 


In this paper, we have attempted to identify the set of chem¬ 
ical reactions that it is essential to include in any chemical 
network used to determine J cr it, the critical UV flux re¬ 
quired to suppress H 2 cooling in atomic cooling halos with 
T v i r > 10 4 K. To do this, we have made use of the reaction- 
based reduction algorithm developed by [Wicbc, Semenov] 
& Henning (2003). This has previously been applied to the 


study of the chemistry of the local interstellar medium, but 
to the best of our knowledge, the present paper marks its 
first use in the study of the chemistry of primordial gas. 

Using this reduction technique with a conservative 
choice of e = 10~ 4 for the reaction weight below which we 
do not retain reactions in our reduced network, we find that 
we can reduce our initial 30 species, ~ 400 reaction network 
to a reduced network containing only eight chemical species 
linked by 26 reactions. We have verified that simulations 
carried out using this reduced network predict essentially 
identical values of Jcrit to those carried out using our full 
network. We have also explored the effect of varying e, and 
find that we continue to be able to predict J cr i t accurately 
as long as e ^ 10 -2 . Setting t = 10 -2 allows us to produce 
an even smaller “minimal” reduced network, containing the 
same set of eight chemical species, but now linked by only 
18 reactions. 

Most of the reactions included in our reduced networks 
are familiar from previous studies of the chemistry of pri¬ 
mordial gas. The main exception is the reaction 

H + H ->• H + + e~ +H. (12) 

This was included in the original investigation by |Omukai] 
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1 2001) of the effect of strong UV radiation fields on the grav¬ 


itational collapse of primordial gas, but has been omitted in 
most subsequent studies. We have investigated the influence 
of this reaction and show that omitting it leads to a 20-30% 
reduction in J cr it- We also confirm the previous finding of 


Latif et al. (2014) regarding the importance of accounting 


for dissociative tunneling when treating H 2 collisional dis¬ 
sociation, and show that omitting this process leads to one 
overestimating J cr i t by almost a factor of two. 

We have compared the values for J cr it that we recover 
using our reduced network with those that we obtain us¬ 
ing several other chemical networks that have previously 
been adopted in studies of the direct collapse model: the re- 
act-xrays network from the krome astrochemistry package, 
the enzo primordial chemistry network, and the |Omukai 1 
1 2001 1 network. In carrying out this comparison we have up¬ 


dated the chemical rate coefficients used in these networks 
to match those used in the current paper, to allow us to 
focus on the uncertainties introduced by differences in the 
set of reactions chosen to construct the different chemical 


networks. We find that the Omukai (2001) network predicts 


Jcrit accurately for the runs with a T4 spectrum, but signif¬ 
icantly overestimates it for runs with a T5 spectrum, owing 
primarily to its neglect of the reaction 


H~+H—i-H + H + e - 


(13) 


We also show that the enzo network tends to predict values 
of Jcrit that are 20-50% larger than those predicted by our 
reduced model, while the KROME react-xrays network pre¬ 
dicts values that are 20-50% smaller. The total uncertainty 
introduced into estimates of J cr it given in the literature due 
to differences between the chemical networks adopted by 
different studies can therefore approach a factor of three. 

To put this number into context, we note that in the 
regime relevant for the direct collapse model, the probability 
of an atomic cooling halo being exposed to a local radiation 
field with a strength J21 > Jcrit is a strongly decreasing func¬ 
tion of Jcrit flDijkstra et al.|2014 l. For example, Inayoshi & 


Tanaka (2014 1 show that for lti 3 < J cr it < 10 4 , this prob 
ability scales approximately as P oc J“ ; ®. A factor of three 
uncertainty in J cr it can therefore potentially correspond to 
a factor of ~ 200 uncertainty in the cosmological number 
density of massive black holes formed by direct collapse. 

Finally, we stress that this uncertainty can be elimi¬ 
nated from future studies of the direct collapse model simply 
by ensuring that the chemical network adopted includes the 
full set of chemical reactions that are important for deter¬ 
mining Jcrit- We therefore recommend that in future work, 
researchers take care to include, at the very least, all of the 
reactions making up the minimal reduced model described 
in Section 1331 
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APPENDIX A: REVISIONS TO THE 
THERMAL MODEL 

We have improved on the thermal model introduced in 
|Glover fe Abel| ( |2008| l and |Glover fe Savin| ( |2009| ) by up¬ 
dating the rates of two of the cooling processes in an effort 
to use the most up-to-date data available. Details of our 
changes are given below. 


A1 Updated cooling rates 

A 1.1 Collisional excitation of H 2 by electrons 


The cooling rate due to collisions between H 2 molecules and 
free electrons that we used in Glover & Savin (20091 was 
taken from Glover & Abel (20081 and was based on rather 
old data from Draine, Roberge & Dalgarno (1983). We have 
now updated our treatment of this process and use a rate 
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based on data from the recent compilation of |Yoon et al.| 
12008). The resulting cooling rate is well fit by the expression 


log Ah 2 ,e = - 21.928796+ 16.815730 log(T 3 ) 

+ 96.743155 log(T 3) 2 
+ 343.19180 log(T 3 ) 3 
+ 734.71651 log(T 3 ) 4 
+ 983.67576 log(T 3 ) 5 
+ 801.812471og(T 3 ) 6 
+ 364.14446 log(T 3 ) 7 
+ 70.609154 log(T 3 ) 8 , 

at temperatures 100 < T < 500 K, and by 

logA h 2 ,e = - 22.921189+ 1.6802758log(T 3 ) 

+ 0.93310622 log(T 3 ) 2 
+ 4.0406627 log(T 3 ) 3 

- 4.7274036 log(T 3 ) 4 

- 8.8077017log(T 3 ) 5 
+ 8.9167183 log(T 3 ) 6 
+ 6.4380698 log(T 3 ) 7 

- 6.3701156 log(T 3 ) 8 , 


(Al) 


(A2) 


at temperatures 500 < T < 10 4 K, where T 3 = T/1000 K. 
These fits are accurate to within 5% over the quoted temper¬ 
ature range. This revised treatment yields less cooling at all 


temperatures than the rate given in Glover & Abel (20081, 


with the effect being particularly pronounced at tempera¬ 
tures around 1000 K. However, the effect on the total H 2 
cooling rate is less significant, since H 2 -H + collisions lead 
to substantially more cooling than H 2 -e” collisions in gas 
where n H + ~ n e . 


A1.2 Collisional excitation of H 2 by protons 

We have also updated our treatment of the cooling rate 
due to collisions between H 2 molecules and protons. The 


treatment in Glover & Savin (20091 was based on Glover & 
Abel| ( 2008|) a nd made use of data from |Gerlich| ( |1990| ) and 
Krstic (20021. Our new treatment makes use of the excita¬ 


tion rates recently calculated by Honvault et al. (2011 2012) 


for the transitions for which these are available, supplement¬ 


ing them with data from Gcrlich (1990) and Krstic (2002) 


for those transitions for which newer data is not available. 
The resulting cooling rate is well fit by the expression 

logA H2iH+ = -22.089523+ 1.5714711 log(T 3 ) 

+ 0.015391166 log(T 3 ) 2 

- 0.23619985 log(T 3 ) 3 

- 0.51002221 logCZs ) 4 

+ 0.32168730 log(T 3 ) 5 , (A3) 

for temperatures in the range 10 < T < 10 4 K. Again, the 
revised treatment yields less cooling at all temperatures than 


the rate given in Glover & Abel (2008). In conditions where 


H 2 -H 4 " collisions make the dominant contribution to the H 2 
cooling function, the total H 2 cooling rate can be as much 
as a factor of two smaller. 


APPENDIX B: REVISIONS TO THE 
CHEMICAL MODEL 

We have improved on the chemical model introduced in 


Glover & Savin (2009) by including three new reactions and 


updating the reaction rates used for a further eight reac¬ 
tions. Details of our changes are given below. 


B1 New reactions 

Bl.l Collisional ionization of atomic hydrogen 

Because the rate at which electrons bring about the colli¬ 
sional ionization of atomic hydrogen 

H + + e“ + e" 


H + e~ 


(Bl) 


is substantially faster than the ionization rate due to H-H 
collisions 


H + H -+ H+ + e" + H, 


(B2) 


the latter process has been omitted from most chemical 
models of primordial gas. However, it can actually dominate 
if the fractional ionization of the gas falls below x ~ 10 —4 , 
or if there is a substantial population of H atoms in excited 
electronic states ( [Om ukai 2001). We therefore include this 
process in our revised chemical model, using the following 
rate coefficient (Lenzuni, Chernoff & Salpeter| 1991; Omukai 

|200l| ). 


fcci.H = 1.2 x 10 _ 17 T 1 ' 2 exp 


' —157800^ 3 _! 

- j, -j cm s . (B3) 


This rate coefficient is based on the experimental cross- 


sections measured by Gealy & van Zyl (19871, and assumes 


that both hydrogen atoms are in their ground state. As we 
discuss in more detail in Paper II, a number of other ver¬ 
sions of the rate coefficient for this reaction can be found 
in the literature (|Drawin||1969| |Hollenbach fe McKee||1979[ 


1989 Kune fe Soon|1991 Soon 1992 Barklem|2007 1 . These 

expressions differ by large amounts at the temperatures rel¬ 
evant for this study and it is unclear which expression gives 
the best description of the actual behaviour of this reaction. 
In this paper, we have chosen to use the version of the rate 
coefficient given in [Lenzuni, Ch ernoff fe Salpeterj (199l[) for 
consistency with the earlier study of Omukai (2001 1 , but in 


Paper II we explore the sensitivity of J cr it to this choice. 

We also include the collisional ionization of atomic hy¬ 
drogen by collisions with neutral helium 


H + He -+ H + + e“ + He. 


(B4) 


We use the rate coefficient computed by |Lenzuni, Chernoff 
& Salpeter ( 1991| ) using cross-sections from Van Zyl, Le & 
Ammel (119811): 


^-cijHe = 1.75 x 10 _ 17 T 1,3 exp 


' — 157800^ 3 _r 

--- J cm s . (B5) 


Again, this expression assumes that all of the hydrogen 
atoms are in their ground state. 

In practice, trapping of Lyman-a photons in the col¬ 
lapsing gas owing to the high optical depth in the Lyman-a 
line leads to the existence of a non-zero population of hydro¬ 
gen atoms in excited electronic states (see e.g. Qmukai |2001 


Schleicher, Spaans & Glover 20101. However, at the densities 


of interest in this study, only a very small fraction (of order 
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10 10 or less) of the total number of hydrogen atoms are 
found in these states (Schleicher, Spaans & Glover 2010), 


and so there is no need to account in the chemical model for 
the effects of collisional ionization out of these states. 


B1.2 Collisional dissociation of H% by atomic hydrogen 

We now account for the destruction of Hi/ by the reaction 

Hj + H -+ H+ + H + H. (B6) 

The rate at which this reaction proceeds depends on the in¬ 
ternal excitation of the Hj molecular ion. At low densities, 
it is safe to assume that all of the ions are in the vibrational 
ground state. In this limit, we adopt the following rate co¬ 
efficient for this reaction 

fccd,o = 1-54 x 10 T exp I———I cm s , (B7) 

which we have derived based on the cross-sections presented 
in 


Krstic & Janev (2003). At high densities, on the other 
hand, the vibrational level populations of the ion ap¬ 
proach their local thermodynamic equilibrium (LTE) values. 
In this limit, we use the LTE rate for this reaction derived by 


Coppola et al. (20111, which was again based on the Krstic 
& Janev (2003) cross-sections. At intermediate densities, we 
interpolate between these two limiting cases by adopting a 
rate coefficient of the form 


/b C d — kc 


k c 




(B8) 


where fc c d,LTE is the rate coefficient in the LTE limit, a = 
(1 + n/ricrit) -1 , and n cr it is the critical density for Hj (i.e. 
the density at which the effects of collisional de-excitation 
and radiative de-excitation become comparable). 

As explained in G lover fe Savin| ( [200 9j), in primordial gas 
the dominant contributions to the collisional excitation or 
de-excitation of H(/ come from collisions with atomic hydro¬ 
gen or with electrons. When collisions with H atoms domi¬ 
nate, we can estimate the critical density by taking the ratio 
of the cooling rates due to Hj in the LTE and low density 
limits, using the expressions for these given in | Glover <%z| 
Savin (2009). The resulting values are reasonably well ap¬ 


proximated by the expression 

ncrit, h — 400r 4 _1 cm -3 , (B9) 

where T 4 = T/IO 4 K. A similar procedure applied to col¬ 
lisions with electrons yields a critical density of electrons 
that does not vary significantly with temperature, and that 
is roughly an order of magnitude smaller than n cr it,H, 


ncrit ,e — 50 cm 


(B10) 


We combine the contributions from atomic hydrogen and 
free electrons by taking a weighted harmonic average, yield¬ 
ing 


^crit 


Xjj 

^-crit.H 


+ 


^crit.e 


(BH) 


Although this is a somewhat crude approximation, it is 
sufficient for our purposes, since J cr it is not very sensi tive to 
the treatment adopted for H)/ excitation (see Section 3.41. 


B2 Updated reaction rates 


Since the publication of Glover & Savin (20091, new exper 


imental and/or theoretical data has become available for a 
number of different reactions. We have therefore updated 
the rates adopted for several of the included reactions, as 
detailed below. 


B2.1 Associative detachment of H with H 
The rate coefficient for the reaction 


IU + H -+ H 2 + e“ 


(B12) 


was recently measured by Kreckel et al. (2010) for tempera¬ 
tures in the range 1 < T < 10 4 K using a merged-beam ap¬ 


proach (see also Bruhns et al.|2010 Miller et al.|20lf I. The 
estimated systematic error in their measurements is around 
25%. Their results disagree with the flowing afterglow re¬ 


sults of Martinez et al. (2009) at 300 K by a factor of two, 


for reasons which remain uncertain, but agree well with the 
measurements made by |Gerlich et al. (2012]) at temperatures 
10 < T < 150 K using an ion trap. We therefore adopt the 
rate coefficient of Kreckel 'et al.| (2010) for this reaction. We 
also adopt the same rate coefficients for the isotopic variants 
of this reaction (i.e. reactions where one or both H atoms 


are replaced by D atoms), following Miller et al. (2012 I, who 


found that there is no significant isotope effect in the cross- 
section for this reaction. 


B2.2 Mutual neutralization of H with H + 

Until relatively recently, the rate of the mutual neutraliza¬ 
tion reaction 


LU + tn 


H + H 


(B13) 


at low temperatures was quite unclear. Glov er, Savin fc| 
Jappsen (20061 surveyed the range of rates given for this re¬ 


action in the astrophysical literature as of 2006, and showed 
that there was almost an order of magnitude scatter in the 
values. Fortunately, the last few years have seen a significant 
improvement in this area. New theo retical (|Stenrup, Larson | 
& Elandcr 2009) and experimental ( Urbain et al.|20~12 l de¬ 
terminations of the rate coefficient yield values that are in 
good agreement with the ones derived by |Croft, Dickinson 


(1986). We therefore adopt the Croft, Dickinson & Gadea 


& Gadea (1999) from the cross-sections of Fussen & Kubach 


(19991 rate coefficient for use in our study. 


B2.3 Dissociative recombination of 


The rate coefficient used in Glover & Savin ([2009) for the 
reaction 


HJ + e~ 


H + H 


(B14) 


was a fit made by|Abel et al. (19971 to the data of Schneider 


et al. (1994 1997|. It assumes that the Hj molecular ions 


are in their vibrational ground state. However, this is true 
only at low densities. At high densities, the H/ - level popu¬ 
lations tend towards their local thermodynamic equilibrium 
(LTE) values, and the resulting dissociative recombination 
rate can be significantly larger (see e.g. the comparison of 
low density and LTE rates in Figure 9 of |Coppola et al.| 
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2011). To account for this, we use a density dependent rate 


coefficient of the form 


k(] r — k. 


dr,LTE 


h 


dr,0 


(B15) 


where fcdr.o is the rate coefficient in the low density limit 
(taken as before from Abel et al.||1997l, fcdr.LTE is the rate 


coefficient in the LT E limit (taken from Coppola et al.|2011 
based on data from Takagi|[2002 1 , a = (1 + n/n CI it) , and 


ricrit is the critical density for Hi, calculated as outlined in 
Section IB1.2I above. 


B2-4 Collisional dissociation of H 2 by H 


We now use the fitting formulae given in |Martin, Schwarz| 
|fc Ma ndy ( 1996]) to determine the rate coefficient for the 
reaction 


H 2 + H-)-H + H + H, 


(B16) 


in place of the combination of data from several sources used 
in Glover fe Savin| ([2009). Importantly, we include the con¬ 
tribution to the total H 2 dissociation rate due to excitation 
to a quasi-bound state followed by dissociative tunneling to 
an unbound state (y d t in the notation of |Martin, Schwarz| 
& Mandy 19961. At low temperatures (T < 4500 K in the 


low density limit), this effect is more important than direct 
collisional dissociation and it is important to account for it 


if one wants to determine J cr it accurately (see Latif et al. 
2014 or Section 3.5 of the present paper). 


B2.5 Charge transfer from H + to H 2 


Glover & Savin (2009) adopted a rate coefficient for the 


reaction 


H 2 + H"* 


H++H 


(B17) 


that was taken from Savin et al. ( 2004a|b I and that assumes 
that all of the H 2 molecules are in their vibrational ground 
state. However, in the LTE limit, the actual rate can be as 
much as an order of magnitude larger (Coppola et al. |201 1). 
In the present study, we therefore adopt a density-dependent 
rate coefficient of the form 


kct — kc 


fc< 


'ct,0 


ct,LTE 


(Bis) 


where fc c t,o is the rate coefficient in the low density limit 
(taken from |Savin et aL||2004a|b| ), fc c t,LTE is the rate coef¬ 
ficient in the LTE limit (taken from Coppola et al. 20111, 
a = (1 + n/ricrit) - , and n cr it is the H 2 critical density. We 
conservatively assume that the latter has the same value 
here as for collisions between H 2 and H. 


B2.6 Three-body formation of H 2 

A large variety of different rate coefficients have been given 
in the astrophysical literature for the reaction 


H + H + H->H 2 + H, 


(B19) 


as summarized in Glover (20081 and Turk et al. (2011 1 . 


Typically, these values were derived by taking the 


experimentally-measured rate for the inverse process (col¬ 
lisional dissociation of H 2 by H) and then applying the prin¬ 
ciple of detailed balance. In general, the resulting rate co¬ 
efficients agree fairly well at high temperatures (where the 
collisional dissociation rate can be easily measured), but dis¬ 
agree substantially at low temperatures, owing in large part 
to the errors introduced by extrapolating the collisional dis¬ 
sociation rates outside of the measured temperature range. 

This unsatisfactory state of affairs was recently ad¬ 
dressed by Forrey (2013a b). He computed a rate coefficient 
for this reaction using a new technique that does not rely 
on detailed balance and that hence should be far more reli¬ 
able at low gas temperatures. The resulting rate coefficient 
is well fit by the simple expression ( Forrey|2013a l 

(B20) 


fc 3b = 6 x 10 _ 32 T _1/4 + 2 x 10~ 31 T “ 1/2 cm 6 s _1 , 
and we use this value in all of our calculations. 


B2 .7 Photodissociation of Hf 
y Glover 

Hj + 7 —» H + + H 


The rate adopted by Glover & Savin (20091 for the process 

(B21) 


assumed at T5 spectrum and also that all of the H^" molec¬ 
ular ions were in their vibrational ground state. The latter 
assumption is reasonable at the redshifts of interest in this 
study provided that the gas density is significantly lower 
than the H^ critical density. However, as n — > n cr it, it be¬ 
comes important to account for the effects of vibrational 
excitation, as this can potentially have a large influence on 
the photodissociation rate (Galli & Palla 1998). 

In our present study, we have therefore calculated H 2 
photodissociation rates for the limiting cases where all of 
the molecular ions are in the v = 0 state and where the 
vibrational level populations have their LTE values. In the 
v = 0 case, we use the expression 


fc pd , H + iV=0 = 2.0 x lO 1 ^ 59 exp (- 


82000 


T ra d ) 


\ s ^ 1 


(B22) 


given by Galli & Palla (19981, based on data from Dunn 


(1968), to compute the photodissociation rate for a black- 


body spectrum and then rescale both the strength of the 
spectrum and the size of the rate by the same amount so 
that the resulting specific intensity at the Lyman limit is 
given by 10 -21 ergs _ 1 cm“ 2 Hz _ 1 sr^ 1 . For a T4 spectrum, 
this procedure yields a photodissociation rate given by 


k 


pd,H+,v =0 ~ 1-74 x 10 J 21 s , 

while for a T5 spectrum we have 


k 


pd,H 


+ n = 5.77 x 10" 

,v=0 


u T — 

J 21 S 


(B23) 


(B24) 


In the LTE case, the Hj photodissociation rate depends on 
the gas temperature, which we cannot assume is the same as 
the radiation temperature. We therefore cannot use the for¬ 


mula given in Galli & Palla (19981, which does assume that 


Tgas = Trad, but instead compute the required rates using 


the cross-sections given in Stancil (1994). These are valid 


for temperatures in the range 3150 < T < 25200 K. We find 
that for temperatures in this range, the H^ photodissocia¬ 
tion rate for a T4 spectrum is well-fit in the temperature 
range 3150 < T < 9000 K by the following expression 
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k 


„ ln _ 8 ( 8500 \ 
= 3.45 x 10 exp-— ) 


J 21 s 


(B25) 


pd,HJ,LTE "■•i' \ j> j 

and at temperatures 9000 < T < 25200 K by the expression 

8500 \ 


k 


pd,H+ LTE 


= 2.5 x 10“ 7 T“°- 22 exp I -- 


T J 


J 21 s~\ (B26) 


Within the range of temperatures quoted above, the fitting 
error of these expressions is no more than 10%. At tempera¬ 
tures T < 3150 K, we extrapolate using Equation [B25| until 


we reach the rate given by Equation |B23| since the v = 0 
rate is also the appropriate low T limit of the LTE rate. 
At temperatures T > 25200 K, we could in principle again 
extrapolate using Equation |B26[ but in practice, the gas in 
our models never reaches these temperatures. 

For the T5 spectrum, we use a similar procedure. In this 
case, the LTE rate is well-fit over the whole temperature 
range 3150 < T < 25200 K by the expression 


k 


pd.H^.LTE 


= 2.2 x 10 -11 T -0 ' 22 J21 s _1 . 


(B27) 


As before, at temperatures T < 3150 K, we extrapolate the 
rate coefficient using the same expression until we reach the 
value given by Equation |B24[ 

Finally, in our numerical models we smoothly interpo¬ 
late between the two limiting cases using the following ex¬ 
pression: 


k 


k 


pd,K 


= k 


pd,H+ v=0 


pd.HJfLTE I £. 


(B28) 


pd,HJ,LTE 


where a = (1 + n/n CT it) 1 , and n CT it is the cr itical density 
of Ht, computed as described in Section 


B1.2 


above. 


fit given in the recent paper by 


Latif et al. (2015 I, and based 


on the work of Coppola et al. (2011) 


fc raH + = dex [-18.20- 3.1941ogT+ 1.786(logT) 2 

-0.2072(logT) 3 ] cm 3 s -1 . (B33) 

The fit gives values for the rate coefficients that agree with 


those tabulated by Ramaker & Peek (1976) to within around 


5-6%. Moreover, at the temperatures most relevant in this 
study, T ~ 4000-8000 K, the agreement is even better, with 
the fit differing from the tabulated values by no more than 
1 - 2 %. 


B2.8 Formation of by radiative association 
Rate coefficients for the reaction 


H + FT 


at +7 


(B29) 


have been computed by both Ramaker & Peek (1976| and 
Stancil, Babb & Dalgarno (19931, and agree to within 3%. 


However, the analytical fits to this data used in most cur¬ 
rent models of primordial chemistry are rather less accu¬ 
rate. Many models use the fit introduced by Shapiro fe Kang| 
(1987), which has the form 


k „+ = 1.85 x io- 23 T 18 cm 3 s _1 


(B30) 


for T < 6700 K and 


k „+ = 5.81 x 10“ 

ra,HT 


(^00) 


V (T) 


3 -1 
cm s 


(B31) 


for T > 6700 K, w here r/(T) = -0 .6 657 lo g(T/56200). On 
the other hand, the Glover &; Savin| p0 09 ) model uses in¬ 
stead the fit given in Galli & Palla ( 1998| ) 

fc ra,H+ = dex [-19.38 - 1.523logT + 1.118(logT) 2 


-0.1269(logT) 3 ] cm 3 s _1 . 


(B32) 


Although both of these analytical fits are based on the [Ra¬ 


maker & Peek (1976) data, they differ from each other, and 


from the rates tabulated in by Ramake r fe Peek) by as much 
as 30%. For this reason, in our current study we do not use 
either of these prescriptions. Instead, we use the analytical 
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